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Quantum chemistry methods exploiting density-functional approximations for short-range 
electron-electron interactions and second-order Mpller-Plesset (MP2) perturbation theory for long- 
range electron-electron interactions have been implemented for periodic systems using Gaussian-type 
basis functions and the local correlation framework. The performance of these range-separated dou¬ 
ble hybrids has been benchmarked on a significant set of systems including rare-gas, molecular, 
ionic, and covalent crystals. The use of spin-component-scaled MP2 for the long-range part has 
been tested as well. The results show that the value of /r = 0.5 bohU 1 for the range-separation 
parameter usually used for molecular systems is also a reasonable choice for solids. Overall, these 
range-separated double hybrids provide a good accuracy for binding energies using basis sets of 
moderate sizes such as cc-pVDZ and aug-cc-pVDZ. 


I. INTRODUCTION 

Among the wide variety of quantum chemistry meth¬ 
ods proposed in the last decade, a great fascination 
resides in the possibility to combine density-functional 
theory (DFT)[1, 2] and explicit many-body correlation 
methods, such as second-order Mpller-Plesset (MP2) 
perturbation theory, [3] random-phase approximations 
(RPA)[4] and coupled-cluster theory. [5] The hope is to 
get the best out of both worlds, that is to combine a 
proper description of long-range dispersion interactions 
(without introduction of empirical corrections [6]) and a 
good description of short-range electron correlations with 
a reduced dependence on the basis set. 

Different ways to combine DFT and wave-function 
techniques have been proposed. These include global 
double-hybrid approaches [7, 8] and range-separated 
approaches. [9-16] In this work we focused on the lat¬ 
ter type, which was initially introduced by Stoll and 
Savin in the 80s. [17] It has been shown that for molecu¬ 
lar complexes (and especially for dimers involving rare- 
gas atoms) it provides excellent performance as regards 
bond lengths, dissociation energies, and harmonic fre¬ 
quencies. In many cases, results are, quote, “superior, 
with medium-size basis sets, to pure DFT and pure 
coupled-cluster calculations” .[14] 

Concerning the study of solids, we note that the in¬ 
troduction of wave-function-based correlation treatment 
for periodic systems is relatively recent. [18-26] Martinez- 
Casado and coworkers proposed to estimate the correla¬ 
tion energy as a MP2 contribution using a B3LYP ref¬ 
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erence state, and applied this scheme to the study of 
the adsorption of helium atoms on a MgO surface. [27] 
Del Ben and coworkers [24] benchmarked the performance 
of some global double hybrids on a set of molecular 
crystals. In particular, these authors used the DSD- 
BLYP functional, [28] that includes a spin-component- 
scaled (SCS)[29] MP2 contribution. Recently, some of us 
applied one-parameter double-hybrid methods to molec¬ 
ular crystals. [30] 

In this work, we tested the performance of range- 
separated double-hybrid methods, [11] combining short- 
range density-functional approximations with long-range 
MP2 correlation, for evaluating the cohesive energy of 
crystalline periodic systems. We implemented this ap¬ 
proach in the Crystal[31] and Cryscor[18] programs, 
that use a basis set of Gaussian-type orbitals centred on 
atoms. The MP2 correlation part is treated within a local 
approach, where the pair-specific virtual space is spanned 
by projected atomic orbitals (PAOs), restricted to the so- 
called domains, i.e. several atoms surrounding the con¬ 
sidered localized occupied orbitals. This technique was 
initially proposed by Pulay in the 80s,[32, 33] and ex¬ 
tended by Werner, Schiitz and coworkers to high-level 
correlated methods for molecules. [34-45] The adaptation 
of the local MP2 method to periodic systems has been 
done in the last decade, [18, 46, 47] and over the years 
has been successfully applied to quite a rich variety of 
systems. [48-60] 

The paper is structured as follows. In Section II, a re¬ 
view of the formal aspects of the methods and their ex¬ 
tension to periodic systems is given. Section III contains 
the details on the model systems, basis sets and compu¬ 
tational parameters used in the calculations. The tests of 
the range-separated double hybrids on three representa- 
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tive systems with the aim of studying the dependence on 
the range-separation parameter are discussed in Section 
IV. Further, Section V presents the benchmarks of sev¬ 
eral range-separated double-hybrid approximations on a 
wider set of systems. Finally, conclusions and perspec¬ 
tives on future work are provided in Section VI. In Ap¬ 
pendix A, the methods employed in this study have been 
additionally tested on molecular dimers cut out from the 
bulk systems. 


II. COMPUTATIONAL METHODS 
A. Periodic range-separated hybrid scheme 

In the range-separated hybrid (RSH) scheme, [11] the 
ground-state energy is approximated with the following 
minimization over (normalized) single-determinant wave 
functions $ 

£rsh = nun {<*|f + V e xt + W*\9>) + -®Hxc N>]} .(!) 

where T is the kinetic energy operator, T4xt is the ex¬ 
ternal potential (nuclei-electron + nuclei-nuclei interac¬ 
tions) operator, FF R is a long-range electron-electron 
interaction operator made with the long-range inter¬ 
action Wg r e (r) = erf(/zr)/r, and -Eg X c[ n 4 > ] is the cor¬ 
responding /z-dependent short-range Hartree-exchange- 
correlation density functional evaluated at the density of 
$. The parameter /z controls the range of the separation. 
It is somewhat clearer to rewrite Eq. (1) as 

Ersh = mini (<b| T + £xt|$> + + *£ hf[^] 

<E> L 

+^XC [ n $] } ! (2) 

where En[n$\ is the usual Hartree energy (with the 
Coulomb electron-electron interaction w ee (r) = 1/r), 
E]JhfW is the long-range Hartree-Fock (HF) ex¬ 
change energy, and K x [.[n$] is the short-range exchange- 
correlation energy. The minimization in Eq. (2) leads 
to familiar hybrid Kohn-Sham (KS)-type self-consistent 
equations determining the RSH orbitals fa and orbital 
energies £j 

F RSH \ fa )= Ei \ fa ), (3) 

with the RSH Fock operator 

pnsu = f + y ext + ^ + ^ + ysr, (4) 

where Vr is the usual Hartree potential operator, V^hf 
is the long-range HF exchange potential operator, and 
V x s c r is the short-range exchange-correlation potential op¬ 
erator. For /z = 0 the RSH scheme reduces to pure KS 
DFT, while for /z —>• oo it reduces to pure HF theory. 

The RSH scheme and other similar range-separated hy¬ 
brid DFT scheme are available in many molecular quan¬ 
tum chemistry programs. A very similar scheme to RSH, 


called RSHX, [61] where the separation is done on the 
exchange energy only, has been implemented for periodic 
systems using a plane-wave/projector-augmented-wave 
(PAW) approach. [62] Another kind of range-separated 
hybrids, called screened-exchange hybrids [63] which use 
short-range HF exchange instead of long-range HF ex¬ 
change, has also been implemented for periodic systems 
using Gaussian-type basis functions [64] or a PAW ap¬ 
proach [65]. Here, we give the main equations for a spin- 
restricted closed-shell RSH scheme for periodic systems 
using local basis functions (see, e.g., Refs. 66-69 for more 
details on periodic HF or KS implementations with lo¬ 
cal basis functions). Due to translational symmetry, the 
crystalline orbitals are labeled by a wave vector k and ex¬ 
panded as \ fa{k)) = E At c M*( k )I^M( k )) where |z/v( k )} = 
N~ 1/2 J2 S e* k ' s |x®) are Bloch functions, and N is the 
number of crystal cells and (r|x®) = X^{ Y ~ g) is an 
atomic-orbital basis function (a contraction of Gaussian- 
type orbital functions) located in the cell characterized 
by the direct lattice vector g. The orbital coefficients 
and energies are found by solving the self-consistent-field 
(SCF) equation at each point k 

F RSH (k)c,(k) = e i (k)S(k)c i (k) l (5) 

with the RSH Fock matrix -F)^ SH (k) = ]T g e ?k ' s .F / RS g lr 
where F^ s g H = (x° |F RSH |xf), and similarly for the over¬ 
lap matrix S IJV (k). The matrix F R ^ & g H is expressed as 


pRSH 
r liv g 


hui'K V J u, 


K] 


V ! 

* \ 


XC,/liZ/g 1 


( 6 ) 


where h^g = (%°|T + Kxt|x®) are the kinetic + external 
potential integrals, V x s c r i/i!/g = (x°|K x s c r |xf) are the short- 
range exchange-correlation potential integrals, is 

the usual Hartree potential contribution calculated with 
Coulombic two-electron integrals 


Vg = E p (X°X^ee| xW +I )> (7) 

Acrml 


and K x * vg is the long-range HF exchange potential contri¬ 
bution calculated with long-range two-electron integrals 

E p ™ (xlxT + >l\x7xl)- (8) 

Acrml 


In these expressions, the density matrix P\ a \ is obtained 
from the occupied orbital coefficients as 

P\ai = - / E cxi(k)cP(k)e(e F - £i (k))e' k - ldk , (9) 

WbzE 


where v is the volume of the Brillouin zone (BZ) and £p 
is the Fermi energy. 

Finally, the RSH energy per unit cell takes the form 


-Ersh = E 


i//jg 


g 


V^g + 2 (A ti/ g 


K 


z^g ) 


10) 
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where the short-range exchange-correlation energy per 
unit cell is, e.g. for generalized-gradient approximations, 

= [ n(r)4c(n(r), Vn(r)) dr, (11) 

J unit cell 

where the integration is over one unit cell and n(r) = 
x“ +1 ( r )x™( r ) is the electron density. In 
this work, we use either the short-range local-density- 
approximation (LDA) exchange-correlation functional 
of Ref. 70 or the short-range Perdew-Burke-Ernzerhof 
(PBE) exchange-correlation functional of Ref. 14 (which 
is a modified version of the one of Ref. 71). The method 
will thus be referred to as RSHLDA or RSHPBE, respec¬ 
tively. 


B. Periodic long-range local second-order 
Mpller-Plesset correction 

The RSH scheme does not contain long-range correla¬ 
tion, but it can be used as a reference for a nonlinear 
Rayleigh-Schrodinger perturbation theory [11, 72, 73] to 
calculate the long-range correlation energy. At second or¬ 
der, the long-range correlation energy is rigorously given 
by a standard MP2 expression evaluated with RSH or¬ 
bitals and orbital energies, and long-range two-electron 
integrals [11, 12]. Here, we give the main equations of the 
long-range local MP2 correction for periodic systems. 

After the periodic RSH calculation, the crystalline 
RSH canonical occupied orbitals are transformed into lo¬ 
calized [74] symmetry-adapted [75] mutually orthogonal 
Wannier functions (WFs). As regards the virtual orbital 
space, mutually nonorthogonal PAOs are constructed by 
projecting the individual atomic-orbital basis functions 
on the virtual space [33]. The long-range first-order 
double-excitation amplitudes XA j b are then obtained by 
iteratively solving the following system of linear equa¬ 
tions [33, 46, 76] 

/A'lr I \ ' J z^RSH rp\r n \ Q 'ylr pRSH l 

lv ia,jb “I" | Pac ^ic,jd °db “I" ^ac -Mc,jd r db f 

(cd)e[ij] l J 

E C \ A Z^RSH rr ilr C 

°ac ^kc,jd °db 

(cd)Gii[j] knearj 

E s -c E r/ tkd A k f H .Shu = o, (12) 

(cd)6u[i] kneari 

where i, j, k refer to WF occupied orbitals, and a, b, 
c, d to PAO virtual orbitals (bold indices combine the 
index within the unit cell and the lattice vector). The 
locality is exploited by restricting the sums over PAO 
pairs (cd) to the pair domain [ij] of PAOs spatially close 
to at least one of the WF i or j, or to the domain it[i] (or 
u[j]) which is the union of all [ik] (or [jk] ) where the sum 
over k is in turn limited to WFs spatially close to i (or 
j). In Eq. (12), AA jb = ( ia |We r e |jb) are the long-range 
two-electron exchange integrals in the WF/PAO basis, 


S'ab is the overlap between PAOs, and pA SH and H 
are elements of the RSH Fock matrix in WF and PAO 
basis, which is obtained by transformation of the Fock 
matrix in the atomic-orbital basis F//1 11 . 

The K\l j b integrals are efficiently evaluated through 
a robust [77] density-fitting scheme, suitably adapted for 
periodic systems. [78-81] By introducing an auxiliary ba¬ 
sis set of Gaussian-type functions - here indicated by 
indices P and Q the integrals are approximated as 

^ia,jb ~ E C( p kee|j b ) + E^^elQ ) d fr 

P Q 

-EC(Pke r e|Q )d% (13) 

PQ 

with fitting coefficients 

<C = E^IV^Q) [ j_ 1 ]q,p > ( 14 ) 

Q 


where [ J 1 ] q p is the Q, P element of the inverse of the 
matrix of Coulomb integrals over the auxiliary functions 
J p _q = (P11/r|Q). In Eqs. (13) and (14), the summa¬ 
tion over fitting functions P and Q is limited to suitable 
local fitting domains. 

The long-range local MP2 correlation energy per unit 
cell is then given as 

*£mP 2= E E *'iajb(2^jb-^b,ja),(15) 

(ij)eP (ab)e[ij] 


where the first sum is over occupied WF pairs (ij) taken 
from a truncated list P, in which the first WF i is lo¬ 
cated in the reference unit cell and the second WF j is 
restricted within a given distance to the first WF i. The 
method obtained after adding the long-range local MP2 
correlation energy to the RSH energy will be referred to 
as RSHLDA+MP2 or RSHPBE+MP2. Obviously, for 
/i = 0, the method reduces to pure LDA [82] or pure 
PBE [83], while for /r = oo it reduces to pure MP2. We 
note that, since exactly the same range-separated Fock 
operator used in the SCF iterations is adopted for the 
evaluation of the long-range local MP2 correlation en¬ 
ergy, no contribution from single excitations [84] arises. 

We also consider the SCS variant [29] of MP2 


z^lr 

- C/ c,SCS-MP2 


E E <»> 

(ij)eP (ab)e[ij] 


(cos + css) aA jb 


r J~l] 

-css G 


ib.ja 


(16) 


where the opposite-spin (OS) and same-spin (SS) coeffi¬ 
cients, taken from the original work [29], are cos =6/5 
and css = 1/3. For molecular systems, the SCS variant 
of MP2 has been shown to significantly improve the ac¬ 
curacy of MP2 for energy differences and properties. For 
compactness of notation, the method will be referred to 
as RSHLDA+SCS or RSHPBE+SCS. 
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(e) Hydrogen 



(f) Lithium Hydride 



(h) Silicon 


cyanide 



(g) Lithium Fluoride 



(i) Silicon Carbide 


Figure 1: Pictorial representation of unit cells of 
crystals used as a benchmark in this work. 


III. COMPUTATIONAL DETAILS 

A. Test systems 

A small but representative set of crystalline systems 
was chosen in order to cover the diverse types of chemi¬ 
cal bond typical of solids (see Figure 1): rare-gas crystals 
(Ne, Ar), molecular solids (CO 2 , HCN, NH 3 ), ionic crys¬ 
tals (LiH, LiF) and covalent semiconductors (Si, SiC). 
Metals were not considered in this study since the per¬ 
turbative nature of the MP2 approach does not allow 
calculations on systems with a very small or zero band 
gap. [85] 

The LiH crystal, being the simplest 3D crystal, has 
been the subject of much attention recently,[19, 86-92] 
and convergence of the total MP2 energy with basis set 
was carefully investigated by some of us. [93] The three 
molecular crystals have been chosen in order to evaluate 
different nature of interactions which contribute to the 
total cohesion energy. This includes dispersive and elec¬ 
trostatic quadrupole-quadrupole interactions[94] in CO 2 


crystal and hydrogen bonds in NH 3 crystal. The HCN 
crystal, instead, represents a simple example of a molec¬ 
ular crystal in which both dispersion and hydrogen bond¬ 
ing are present. These particular systems have been stud¬ 
ied by different authors with different approaches, includ¬ 
ing periodic post-Hartree-Fock methods in the last few 
years. [92, 95-100] Rare-gas crystals are of general interest 
as purely dispersion-bonded crystals, where the dominant 
role of electron correlation effects is well known. [101 105] 
Details on the systems and the geometries we used in 
the calculations are reported in Table I. The experimen¬ 
tal lattice parameters, as indicated, were adopted in all 
cases. For molecular crystals, internal coordinates were 
taken from Ref. 95, where they were obtained by a fixed- 
volume optimization of internal coordinates performed at 
the B3LYP-D* level.[117] 


B. Parameters of the calculations 

All the periodic calculations were performed with 
development versions of the Crystal14[31] and 
Cryscor[18, 46, 47, 84] programs. Molecular calcula¬ 
tions were performed either with the above codes or with 
Molpro.[118, 119] As for the Crystal calculations, we 
adopted a homogeneous 8 x 8 x 8 k-point sampling of the 
reciprocal space, and integral-screening tolerances set to 
10 _8 ,10 _8 ,10 _8 ,10 _2O ,10 _50 . For the meaning of these 
thresholds we address the user to the Crystal14 user’s 
manual; [ 120 ] here we just point out that, as in our pre¬ 
vious works, [93, 95] we tightened the thresholds for the 
exchange integrals (last two numbers) with respect to 
defaults. 

In CRYSCOR, the fundamental input parameters refer 
to the locality ansatz. In particular the most relevant one 
is the selection of excitation domains assigned to each oc¬ 
cupied orbital, since in the present calculations the re¬ 
cently implemented orbital-specific-virtual method [ 121 ] 
was not adopted. Domains for rare-gas and ionic solids 
were chosen in order to consider PAOs belonging to two 
coordinated shells of atoms around the occupied orbital, 
while for molecular crystals the domains coincide with 
the molecule to which the orbital belongs. For semicon¬ 
ductor crystals, we selected the domains in order to take 
into account six tetrahedral units for a total amount of 
26 atoms around each bond orbital. The maximum pair 
distance considered for the evaluation of MP2 correlation 
energy was fixed at 12 A. Two-electron integrals within 
this range are evaluated efficiently in CRYSCOR either by 
density fitting [78, 80, 81] or multipolar expansion tech¬ 
niques - if the inter-orbital distance exceeds 8 A. 


C. Basis sets 

As mentioned in the introduction, the Crystal and 
CRYSCOR programs use a basis set of Gaussian-type or¬ 
bital functions centred on atoms to create atomic or- 
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Table I: Structural information about the crystals used in this work. n a to is the total number of atoms per cell. In 
last column we report reference to the works each structure was taken from. 



System 

a; b\ c (A) 

Space group 

^ato 

Ref. 

Rare-gas 

Ne 

Ar 

4.464 

5.300 

Fm3m 

Fm3m 

1 

1 

[106-108] 
[109, 110] 


C0 2 

5.54 

Pa3 

12 

[111] 

Molecular 

nh 3 

5.048 

P2i3 

16 

[112] 


HCN 

4.13; 4.85; 4.34 

Imm2 

3 

[113] 

Ionic 

LiH 

LiF 

4.084 

4.010 

Fm3m 

Fm3m 

2 

2 

[89] 

[114] 

Semiconductor 

Si 

5.430 

Fd3m 

2 

[115] 

SiCJ 

4.358 

F43m 

2 

[116] 


bitals. In this study we mainly employed Dunning’s cc- 
pVXZ.[122] Where possible, unmodified basis sets with 
X=D,T,Q, have been used. 

Exceptions are: [123] 

- LiH and LiF crystals: a suitably optimized cc- 
pVQZ basis set (See supplementary information for 
details) was adopted for Li in all calculations, re¬ 
gardless of the basis on H or F. The basis sets in¬ 
dicated in the tables and figures refer, therefore, 
only to the latter atomic species, for which we used 
standard Dunning’s basis sets. 

- Si and SiC crystals: the cc-pVDZ basis for Si 
needed to be re-optimized in order to allow for SCF 
convergence. However, it was possible to use an un¬ 
modified cc-pVDZ basis on C. 

Augmentation was made only for polarization shells, i.e. 
no diffuse s-type functions were included in any case. 
Hereafter, the prefix p-aug will be adopted. This aug¬ 
mentation scheme proved, in a number of cases, [30, 95- 
97] to be very effective in keeping the beneficial effects 
of polarization-augmented basis sets on the correlation 
energy while moderating the impact on the computa¬ 
tional demands. No dual basis-set scheme[84] is adopted 
in the present work, which would imply contributions 
from single excitations in the MP2 treatment[30]. As a 
consequence, the single-excitation contribution is always 
zero in the calculations presented here. We did not at¬ 
tempt to extrapolate our results to the complete-basis-set 
limit (a recent work [124] suggests an exponential con¬ 
vergence of the long-range correlation energy upon the 
cardinal number of the basis) since this would go beyond 
the scopes of the present work. 

D. Cohesive energies 

We computed the cohesive energy E coh (V) per X 
(where X can stand for either a molecule or an atom) 


at a given volume V of the unit cell, as: 

E coh = ^-E^r ] , ( 17 ) 

where Z is either the number of molecular units or the 
number of atoms in the unit cell; ifbuik and E are 
respectively the total energy per unit cell of the bulk 
system and the total energy of X in the gas phase geom¬ 
etry - which has been optimized at the B3LYP-D* level 
for molecular systems, for consistency with the adopted 
crystalline structures. When X refers to a single atom, as 
in the case of rare-gas crystals, obviously E . 

In order to correct for the basis set superposition error 
(BSSE) we adopted the standard Boys-Bernardi counter¬ 
poise (CP) method: [125] 

^ c C oh = ^coh + ^ ulk] - £^5, (18) 

where F?j£ ulk ] and are the energies of X in the 

crystalline bulk geometry without and with ghost func¬ 
tions, respectively. Note that in the case of ionic crystals, 
cohesive energies are evaluated with respect to the iso¬ 
lated atoms, not ions. The energies for the latter were 
computed in the framework of the spin-unrestricted for¬ 
malism, using the Molpro code. [118, 119] 

IV. DEPENDENCE ON THE 
RANGE-SEPARATION PARAMETER 

The literature discussing the implementations of range- 
separated DFT approaches for molecular systems often 
adopts a value of /i = 0.5 bohr -1 (if not otherwise spec¬ 
ified, bohr -1 will be dropped in the following sections) 
for the range-separation parameter (see Section II). This 
originates from some benchmark on molecular systems, 
and is justified under the consideration that fj, should cor¬ 
respond to the inverse of the average distance between 
valence electrons (i.e. twice the Seitz radius, giving 
around 1-2 bohr in valence regions).[61] 
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dz O qz — B — a-tz — a— exp- 

tz —A— a-dz — • — a-qz —■— 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

|x (bohr' 1 ) 



|i (bohr' 1 ) 


Figure 2: Cohesive energy of the Ne, CO 2 , and LiH 
crystals calculated with the RSHPBE+MP2 method as 
a function of the range-separation parameter fj, for 
different basis sets. 


dz — 0— tz — A— a-dz — •— a-tz — a— exp- 



p (bohr' 1 ) 


Figure 3: Cohesive energy of the CO 2 crystal calculated 
with the RSHPBE+MP2 method as a function of the 
range-separation parameter /1 for different basis sets. 


Solids exhibit a wider variety of chemical bonds with 
respect to molecular complexes: dispersion interactions, 
ionic, covalent, and metallic bond (the latter we can not 
deal with, using the approaches presented in this work). 
The crystalline environment is, by itself, different in na¬ 
ture from that of a molecular system, due to its infinite 
character and close packing of atoms. For these reasons, 
we studied the dependence on the value of /r. The range 
we considered was 0.1-1.0 (corresponding to distances 
between electrons ranging from about 5 to 0.5 A) which 
represents a range of physically reasonable distances. 

In Figure 2 the cohesive energy calculated with the 
RSHPBE+MP2 method is reported as a function of fi 
for three systems of our benchmark set: the Ne, CO 2 , 


Figure 4: Decomposition of the RSHPBE+MP2 
cohesive energy of the Ne, CO 2 , and LiH crystals in 
individual RSHPBE and MP2 contributions as a 
function of the range-separation parameter /i for the 
p-aug-cc-pVDZ basis set. 


and LiH crystals. For each system, curves obtained with 
different basis sets are reported, as well as the experi¬ 
mental value. 

In the case of the purely dispersion-bonded Ne crys¬ 
tal, the cohesive energy is quite dependent on the basis 
set. It stems from the expected strong basis set depen¬ 
dence of the MP2 part, but also from the basis set de¬ 
pendence of the DFT part which is even stronger for this 
system. The latter can be probably explained by the 
very small scale of the interaction energy (of an order 
of just a kJ/mol) and a delicate balance between dis¬ 
persion attraction and exchange repulsion which needs 
to be captured in a proper way. At the same time, for 
any given basis set, the energy is almost independent of 
/r above p, > 0.5. In particular, with the largest basis 
set (p-aug-cc-pVQZ), RSHPBE+MP2 gives an accurate 
cohesive energy for fj, > 0.4. 

For the other two systems, the basis set dependence 
is virtually negligible at the scale of the interaction en¬ 
ergy, provided the polarization-augmented basis sets are 
used. For the CO 2 crystal (middle panel of Figure 2), 
the curves of the cohesive energy exhibit a minimum at 
H ss 0.4 — 0.5 for all basis sets. Near this minimum, 
RSHPBE+MP2 slightly overbinds with the augmented 
basis sets. The obtained cohesive energy is much more 
accurate than pure PBE (/r = 0), as shown in Figure 3 
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where we reported the same curves in an extended range 
of /r. The effect of the diffuse basis functions on the RSH- 
PBE+MP2 cohesive energy is important for /r > 0.4 and 
seems converged with the p-aug-cc-pVDZ basis set. For 
jx < 0.2 the role of many-body correlation effects, not 
included in the exchange-correlation functional, becomes 
insignificant and the impact of augmentation of the basis 
set becomes very small. For the LiH crystal, the cohe¬ 
sive energy curves are quite flat between /r = 0.1 and 0.5, 
with a minimum at around fx « 0.2. Near the minimum, 
RSHPBE+MP2 gives an accurate cohesive energy with 
basis sets larger than the cc-pVDZ basis set. 

In Figure 4 the decomposition of the 
RSHPBE+MP2/p-aug-cc-pVDZ cohesive energy curves 
of Figure 2 in individual RSHPBE and MP2 contribu¬ 
tions is reported. It is interesting to observe how the 
absolute value of the MP2 contribution increases for Ne 
and CO 2 when increasing /i from 0.1 to about 0.5, and 
then saturates. For LiH, the absolute value of the MP2 
contribution increases almost linearly as a function of 
H. This can be rationalized in the following way. By 
increasing the /1 parameter, one progressively includes 
more short-range MP2 correlation. For small values 
of fj, it still means adding more dispersion, leading to 
progressive growth of the magnitude of the attractive 
MP2 contribution. However, at some point, /1 becomes 
so large, that the short-range intra-molecular MP2 
component of the interaction energy, which is usually 
repulsive in molecular crystals, [99] starts contributing. 
Further on it even outruns the short-range dispersion, 
whose accumulation slows down due to the packing 
effects, and the MP2 contribution curve for Ne and CO 2 
turns into a slightly decaying regime. In contrast to the 
molecular crystals, for LiH the short-range correlation 
is stabilizing (correlation substantially strengthens 
binding even in LiH molecule [89]), so in this system the 
magnitude of the MP2 contribution always grows with 
increase of the /r parameter. 

From the present results, we conclude that the value 
of fx = 0.5 is reasonable for the solids considered. Hence, 
as a conservative choice based on previous experience 
on molecular systems, we chose this value of fx for the 
wider benchmark of the range-separated double hybrids 
carried out in Section V. Notably, for that value of /x, 
range-separated double hybrids show a similar basis set 
dependence as full MP2 for the basis sets considered here. 


V. A WIDER BENCHMARK OF 
RANGE-SEPARATED DOUBLE HYBRIDS IN 
SOLIDS 

The results reported in Section IV give a first hint that 
range-separated double hybrids might provide cohesive 
energy of solids that are more accurate than the ones ob¬ 
tained by pure DFT or pure MP2, as generally observed 
for molecular systems. In this Section, we further bench- 
marked the range-separated double hybrids, with double- 


zeta quality basis sets, on the set of solids described in 
Section III A. In Table II cohesive energies are compiled 
for the cc-pVDZ basis set and the following methods: 

- Pure DFT (LDA, PBE), corresponding to fx = 0; 

- Pure MP2 and SCS-MP2, corresponding to /i = 00 ; 

- Range-separated combinations of the above with 
H = 0.5. 

For each system the zero-point-energy corrected experi¬ 
mental cohesive energy is also provided, and a number of 
statistical indicators are calculated. 

Looking at the mean errors (MEs), the mean absolute 
errors (MAEs), or the error variances (er 2 ), it is seen that 
the range-separated double hybrids globally perform bet¬ 
ter for cohesive energies than the respective pure DFT 
and pure MP2 methods. The mean absolute relative er¬ 
rors (MAREs) also confirm this trend, with the exception 
of pure LDA which gives better cohesive energies for Al¬ 
and CO 2 with the cc-pVDZ basis set compared to the 
other methods, leading to the smallest MARE despite 
being less accurate than the other methods for ionic and 
semiconducting crystals. As a matter of fact, the large 
difference in behavior on individual systems observed 
between pure LDA and PBE is greatly smoothened in 
the range-separated double hybrids which yield similar 
results for the two density functionals. Yet, on aver¬ 
age, the PBE-based range-separated double hybrids gives 
slightly more accurate cohesive energies than the LDA- 
based ones. The SCS variants of the range-separated 
double hybrids perform better that the non-SCS ones for 
the semiconducting crystals, but worse for the rare-gas 
and molecular crystals. 

Results for the p-aug-cc-pVDZ basis are analogously 
reported in Table III. Similar trends are observed even 
though the results are more mixed. In comparison to 
pure PBE and pure MP2, the RSHPBE+MP2 method 
provides cohesive energies that are more accurate for Ne, 
Ar, CO 2 , NH 3 , and about as accurate for LiH, LiF, Si 
and SiC. The HCN crystal is a challenging system: it 
is known that pure MP2 (with large basis sets) tends 
to overbind this crystal, [99] and the RSHLDA+MP2 or 
RSHPBE+MP2 method seems to accentuate this behav¬ 
ior. 

Concerning molecular crystals, results can be com¬ 
pared with density-scaled double hybrids reported in 
Ref. 30. For HCN, the comparison is more straightfor¬ 
ward since the same basis set was used in both cases. 
With the DSIDH-PBEsol (A = 0.80) functional, the 
computed cohesive energies are -33.6 and -39.3 kJ/mol 
for the cc-pVDZ and p-aug-cc-pVDZ basis sets, re¬ 
spectively, remarkably smaller than RSHLDA+MP2 and 
RSHPBE+MP2 methods. A similar behaviour is ob¬ 
served for CO 2 and NH 3 , even if the basis sets are not 
exactly the same. For DS1DH-PBE, the cohesive ener¬ 
gies are even smaller. 


Table II: Benchmark of range-separated double-hybrid methods on cohesive energies (-E C oh in kJ/mol) of crystalline 
systems using double-zeta quality basis sets. The RSHLDA+MP2, RSHPBE+MP2, RSHLDA+SCS, and 
RSHPBE+SCS results are obtained with a value of the range-separation parameter oi fi = 0.5. The statistical 
indicators calculated are: mean errors (MEs), error variances (a 2 ), mean absolute error (MAEs), and mean absolute 

relative errors (MAREs). 


Crystal 

LDA 

PBE 

RSHLDA+MP2 

RSHPBE+MP2 

RSHLDA+SCS 

RSHPBE+SCS 

MP2 

SCS „ 

Ref.“) 

^coh 

% 

-^coh 

% 

^coh 

% 

-^coh 

% 

-^coh 

% 

-^coh 

% 

-^coh 

% 

-^coh 

% 

Ne 

-0.31 

84 

-0.27 

86 

+0.11 

106 

+0.10 

105 

+0.18 

109 

+0.18 

109 

+0.22 

111 

+0.28 

114 -1.97 

[107] 

Ar 

00 

cb 

12 

+2.36 

131 

+1.10 

114 

+0.87 

111 

+2.00 

126 

+1.77 

123 

+3.13 

140 

+4.07 

153 -7.73 

[110] 

co 2 

-27.8 

11 

-3.8 

88 

-21.9 

30 

-22.6 

27 

-17.6 

43 

-18.2 

41 

-11.5 

63 

-7.9 

75 -31.1 

[126] 

nh 3 

-55.2 

52 

-26.9 

26 

-32.1 

12 

-32.8 

10 

-29.1 

20 

-29.7 

18 

-24.2 

33 

-19.5 

46 -36.3 

[126] 

HCN 

-51.4 

21 

-28.0 

34 

-40.1 

6 

-40.4 

5 

-36.3 

15 

-36.6 

14 

-31.7 

26 

-27.2 

36 -42.6 

[126] 

LiH 

-265 

10 

-236 

2 

-245 

2 

-238 

1 

-244 

2 

-237 

1 

-226 

6 

-229 

5 -240 

[127] 

LiF 

-482 

12 

-425 

1 

-454 

6 

-438 

2 

-452 

5 

-438 

2 

-432 

0 

-427 

1 -430 

[127] 

Si 

-506 

12 

-433 

4 

-493 

9 

-490 

8 

-469 

4 

-465 

3 

-433 

4 

-404 

11 -452 

[127] 

SiC 

-697 

12 

-603 

4 

-648 

4 

-638 

2 

-633 

0 

-623 

0 

-585 

6 

-559 

11 -625 

[127] 

CT 2 

12.1 


5.0 


6.0 


4.8 


3.9 


2.7 


6.0 


9.9 



ME 

-25.0 


12.6 


-7.4 


-3.6 


-1.3 


2.4 


14.1 


21.9 



MAE 

26.3 


12.6 


13.3 


9.5 


10.0 


7.0 


14.5 


21.9 



MARE 


25 


42 


32 


30 


36 


35 


43 


50 



a1 The experimental sublimation energies of molecular crystals taken from Ref. 126 were corrected for zero-point energy (ZPE) and thermal 
effects at 298 K by a constant 2 RT contribution. [128] Atomization energies of semiconductors and ionic crystals were corrected for ZPE in 
accordance to the zero-point anharmonic expansion correction (values from Ref. 129) 


Table III: Same as Table II with polarization-augmented double-zeta quality basis sets. 


Crystal 

LDA 

PBE 

RSHLDA+MP2 

RSHPBE+MP2 

RSHLDA+SCS 

RSHPBE+SCS 

MP2 

SCS „ 

Ref.“) 

-^coh 

% 

-^coh 

% 

-^coh 

% 

-^coh 

% 

-^coh 

% 

-^coh 

% 

-^coh 

% 

-^coh 

% 

Ne 

-0.44 

78 

-0.40 

86 

-1.24 

37 

-1.23 

38 

-0.87 

56 

-0.86 

56 

-1.10 

44 

-0.74 

62 -1.97 

[107] 

Ar 

-11.9 

54 

+0.42 

105 

-7.50 

3 

-7.65 

1 

-4.88 

37 

-5.03 

35 

-6.45 

17 

-3.57 

54 -7.73 

[110] 

co 2 

-32.9 

6 

-5.8 

81 

-33.6 

8 

-34.4 

11 

-17.6 

43 

-27.5 

12 

-25.2 

19 

-18.5 

41 -31.1 

[126] 

nh 3 

-55.6 

53 

-26.2 

28 

-38.8 

7 

-39.7 

9 

-29.1 

20 

-35.0 

4 

-31.8 

12 

-25.4 

30 -36.3 

[126] 

HCN 

-54.9 

29 

-29.7 

30 

-48.4 

14 

-48.7 

14 

-36.3 

15 

-43.4 

2 

-41.4 

3 

-35.0 

18 -42.6 

[126] 

LiH 

-263 

10 

-233 

3 

-243 

1 

-236 

2 

-241 

0 

-235 

2 

-225 

6 

-229 

5 -240 

[127] 

LiF 

-482 

12 

-426 

1 

-454 

6 

-440 

2 

-453 

5 

-440 

2 

-435 

1 

-429 

0 -430 

[127] 

Si 

-509 

13 

-436 

4 

-474 

5 

-474 

5 

-453 

0 

-448 

1 

-425 

6 

-395 

13 -452 

[127] 

SiC 

-701 

12 

-607 

3 

-652 

4 

-642 

3 

-638 

2 

-627 

0 

-599 

4 

-577 

8 -625 

[127] 

<r 2 

12.6 


4.5 


4.8 


3.4 


3.0 


1.4 


4.6 


8.6 



ME 

-27.1 


11.4 


-9.5 


-6.3 


-3.1 


0.5 


8.5 


17.1 



MAE 

27.5 


11.4 


9.8 


7.4 


5.4 


3.4 


9.6 


17.1 



MARE 


30 


37 


9 


9 


13 


13 


12 


25 



See footnote in Table II. 


VI. CONCLUSIONS 

We implemented and tested range-separated double¬ 
hybrids methods in the Crystal and CRYSCOR pro¬ 
grams for the study of crystalline solids. The approaches 
considered include either LDA- or PBE-type density 
functionals for short-range electron-electron interactions, 
and a local MP2 correlation correction for the long-range 
electron-electron interactions either in its standard or its 
SCS form. 

The value of fi = 0.5 bohr -1 for the range-separation 
parameter commonly adopted for molecular systems was 
found to be also a reasonable choice for solids and was 
thus adopted for this study. The range-separated double¬ 
hybrids methods have been tested on a significant test set 


of cohesive energies of nine prototypical crystalline sys¬ 
tems. A summary of the results is provided in Figure 5. 
With double-zeta correlation consistent basis sets, either 
augmented with diffuse polarization functions or not, the 
range-separated double hybrids are globally more accu¬ 
rate than the respective pure DFT or MP2 calculations. 
As for pure MP2, the effect of augmentation of the basis 
set with diffuse polarization functions is important for 
the range-separated double hybrids, reducing the MARE 
values by about a factor of three with respect to non- 
augmented results. Overall, the SCS variants of the 
range-separated double-hybrids appear to be less accu¬ 
rate than the non-SCS ones, but the reader shall be aware 
that the present benchmark set does not include stacking- 
type interactions, where SCS is expected to bring more 
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Figure 5: Summary of the performance of methods 
tested in this work on crystalline system, measured by 
their MAREs (cf. Tables II and III). Empty lined bars 
refer to the cc-pVDZ basis set and fully colored bars 
refer to the p-aug-cc-pVDZ basis set. 


significant improvements. 

Future developments of range-separated DFT ap¬ 
proaches for solids might include the implementation of 
other flavors of the methods (such as different short-range 
density functionals [130] or different long-range electron 
correlation methods [131, 132]). In particular, the notori¬ 
ous overestimation of dispersion by MP2 in highly polar¬ 
izable systems, which can affect, as observed in this work, 
also the range-separated double-hybrids employing MP2 
as the long-range model, can be cured by substituting the 
latter with approximate coupled-cluster models, contain¬ 
ing only the low-order slowly decaying terms. [133-135] 

We plan also to gain more experience on calculations 
tackling “real-life” problems. The performance on di¬ 
verse properties other than cohesive energies will be of 
interest. One notable example is the relative stability 
of crystalline polymorphs, [49, 55, 136] which tradition¬ 
ally represent a tough challenge for quantum chemistry 
methods. Another quantity of interest, where the range- 
separated double hybrids can become particularly effec¬ 
tive, is physisorption on surfaces or in porous crystals. 
Indeed, since the short-range part of the interaction, in¬ 
cluding the intra-host and intra-adsorbate components, is 
in this scheme described by DFT, the MP2 treatment can 
be reduced to the inter-host-adsorbate pairs only,[137] 
making such calculations computationally very efficient. 
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Appendix A: Molecular dimers cut from the bulk 

In order to highlight the peculiarities of the crystalline 
case with respect to the molecular one, we have per¬ 
formed the same kind of calculations as presented in Ta¬ 
bles II and III on molecular dimers. Since our aim is 
not to discuss the performance of range-separated double 
hybrids on molecular complexes, which has been widely 
explored in the literature, we have extracted dimers from 
the bulk structures without re-optimizing the geome¬ 
tries. This has been done for rare-gas, ionic, and molec¬ 
ular crystals. As a reference, we have performed calcula¬ 
tions at the same geometry with coupled cluster singles, 
doubles, and perturbative triples (CCSD(T)) using an 
p-aug-cc-pV5Z basis set. 

Tables IV and V compile the results for the cc-pVDZ 
and p-aug-cc-pVDZ basis sets. In the case of HCN 
three different types of dimer were chosen, to highlight 
the different types of interaction present in the cluster. 
In Figure 6 a portion of the bulk crystal is reported, 
where some monomers are labeled: one purely hydrogen- 
bonded (“hb”), and two dispersion-bonded ( “db” and 
“db2” ). In the case of LiH and LiF the cohesive energy 
is computed with respect to the neutral atoms. 

Although the test set is not the same as in Sec¬ 
tion V, and is strongly biased towards molecular and 
dispersion-bonded systems, some comparisons on MAEs 
and MAREs can be made. At a first glance it clearly 
appears that these two indicators are significantly higher 
than in the bulk systems. 

For the rare-gas dimers, the results are tremendously 
improved by the addition of diffuse functions in the basis 
set. This effect is almost exactly parallel to what was ob¬ 
served for bulk crystals, and the same can be said about 
the performance of different methods for these systems. 

The molecular dimers show a completely different pic¬ 
ture: here relative errors are quite large in all cases, so 
that these errors dominate the overall statistics. It is in¬ 
deed clear, from inspection of different dimers from the 
HCN crystal structure, how good (excellent, in the case 
of pure LDA) results for the bulk come from error cancel- 
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Table IV: Molecular complexes: benchmark of range-separated double-hybrid methods on interaction energies (-Bint 
in kJ/mol) of dimers cut out from bulk crystals using the cc-pVDZ basis set. The RSHLDA+MP2, RSHPBE+MP2, 
RSHLDA+SCS, and RSHPBE+SCS results are for a value of the range-separation parameter of p, = 0.5. Reference 
CCSD(T) calculations are for an aug-cc-pV5Z basis set. “hb”, “db” and “db2” label three different HCN dimers, 

cfr. Figure 6 


System 

LDA 

PBE 

RSHLDA+MP2 

RSHPBE+MP2 

RSHLDA+SCS 

RSHPBE+SCS 

MP2 

scs 


CCSD(T) 

E'int 

Vo 

Lint 

Vo 

E'int 

% 

E'int 

Vo 

Lint 

Vo 

Lint 

Vo 

Lint 

Vo 

Lint 

% 

[ Ne J 2 

-0.027 

85 

-0.024 

86 

+0.013 

107 

+0.012 

107 

+0.018 

110 

+0.018 

110 

+0.022 

112 

+0.026 

115 

-0.178 

[ Ar i 2 

-0.64 

19 

+0.06 

112 

+0.15 

128 

+0.13 

124 

+0.21 

139 

+0.19 

136 

+0.32 

159 

+0.38 

171 

-0.54 

[C0 2 ] 2 

-4.7 

116 

-0.9 

59 

-2.9 

30 

-3.0 

35 

-2.2 

2 

-2.4 

7 

-1.4 

34 

-0.9 

59 

-2.2 

[nh 3 ] 2 

-18.3 

219 

-10.7 

87 

-9.9 

73 

-10.0 

74 

-9.4 

63 

-9.4 

64 

-7.8 

36 

-6.8 

18 

-5.7 

[HCN] 2 / hb 

-27.0 

173 

-17.5 

76 

-20.4 

106 

-20.5 

107 

-19.5 

97 

-19.6 

98 

-16.3 

64 

-14.6 

48 

-9.9 

[HCNj 2 / db 

+3.8 

74 

+4.3 

95 

+4.5 

105 

+4.4 

101 

+4.8 

119 

+4.7 

115 

+4.3 

96 

+4.5 

107 

+2.2 

[HCN] 2 / db2 

-0.9 

321 

+2.1 

418 

+ 1.8 

349 

+1.7 

310 

+2.4 

484 

+2.2 

446 

+2.3 

467 

+2.8 

608 

+0.4 

Li-H 

-715 

11 

-712 

11 

-702 

10 

-703 

10 

-702 

10 

-703 

10 

-697 

9 

-698 

9 

-641 

Li-F 

-767 

13 

-762 

13 

-745 

10 

-745 

10 

-745 

10 

-745 

10 

-737 

9 

-737 

9 

-676 


13.2 


12.4 


10.3 


10.4 


10.3 


10.4 


9.3 


9.3 



ME 

-21.8 


-18.5 


-15.8 


-16.0 


-15.5 


-15.7 


-13.7 


-13.3 



MAE 

MARE 

22.2 

115 

19.4 

106 

16.6 

102 

16.7 

98 

16.5 

115 

16.6 

111 

14.6 

110 

14.6 

127 




Table V: Molecular complexes: 

same 

as Table IV with the p-aug-cc-pVDZ basis set. 




System 

LDA 

PBE 

RSHLDA+MP2 

RSHPBE+MP2 

RSHLDA+SCS 

RSHPBE+SCS 

MP2 

SCS 

CCSD(T) 

-Lint 

yr 

Lint 

Vo 

Lint 

Vo 

Lint 

Vo 

Lint 

Vo 

Lint 

Vo 

Lint 

Vo 

Lint 

Vo 

i Ne J 2 

-0.038 

79 

-0.035 

80 

-0.080 

55 

-0.081 

55 

-0.055 

69 

-0.055 

69 

-0.072 

59 

-0.047 

74 

-0.178 

[Arj 2 

-1.07 

98 

-0.14 

74 

-0.34 

37 

-0.37 

32 

-0.17 

69 

-0.19 

64 

-0.22 

59 

-0.03 

95 

-0.54 

[C0 2 ] 2 

-5.8 

164 

-1.4 

38 

-4.6 

110 

-4.8 

116 

-3.6 

66 

-3.8 

72 

-3.4 

56 

-2.4 

11 

-2.2 

[nh 3 ] 2 

-17.9 

212 

-10.2 

78 

-11.1 

94 

-11.2 

95 

-10.3 

79 

-10.3 

80 

-9.4 

64 

-8.0 

40 

-5.7 

[HCN] 2 / hb 

-27.9 

182 

-18.1 

83 

-22.1 

123 

-22.1 

123 

-20.9 

111 

-21.0 

112 

-18.1 

82 

-16.2 

63 

-9.9 

[HCN] 2 / db 

+3.7 

70 

+4.3 

99 

+4.3 

98 

+4.2 

94 

+4.8 

118 

+4.7 

114 

+3.7 

68 

+4.1 

88 

+2.2 

[HCN] 2 / db2 

-1.5 

468 

+ 1.9 

362 

+0.7 

66 

+0.5 

19 

+1.5 

260 

+ 1.3 

215 

+0.8 

100 

+ 1.7 

316 

+0.4 

Li-H 

-720 

12 

-717 

12 

-707 

10 

-708 

10 

-707 

10 

-708 

10 

-700 

9 

-700 

9 

-641 

Li-F 

-765 

13 

-760 

12 

-744 

10 

-745 

10 

-744 

10 

-744 

10 

-736 

9 

-736 

9 

-676 

a 1 2 3 4 5 6 7 8 

13.4 


12.6 


10.6 


10.7 


10.6 


10.7 


9.4 


9.4 



ME 

-22.5 


-19.1 


-16.9 


-17.1 


-16.5 


-16.7 


-14.7 


-14.1 



MAE 

MARE 

22.9 

144 

19.8 

93 

17.4 

67 

17.5 

62 

17.2 

88 

17.3 

83 

14.9 

56 

14.6 

78 



lations - overestimation of hydrogen bonds and underes¬ 
timation of dispersion interactions. One must not be sur¬ 
prised to find positive value for the interaction between 
some dimers, because they are not at the equilibrium ge¬ 
ometry. 

Looking at ionic systems - LiH and LiF - we see that 
the results seem almost insensitive to the method and 


to the augmentation of the basis set. We attribute this 
behavior to the strong, extremely short-ranged character 
of the ionic interaction. At the same time, recalling the 
results from Section V, long-range electron correlation 
effects (dispersion) are key to a correct description of the 
bulk. 
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